/*
 *
 * gengraph - generation of random simple connected graphs with prescribed
 *            degree sequence
 *
 * Copyright (C) 2006  Fabien Viger
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
// Pascalou ...
#ifdef pascalou
#define my_random() random()
#define MY_RAND_MAX 0x7FFFFFFF
#else
#include "gengraph_definitions.h"
#endif

#include "gengraph_powerlaw.h"
#include <cstdio>
#include <cmath>
#include <cassert>

#include "igraph_error.h"

namespace gengraph {

// Destructor
powerlaw::~powerlaw() {
  delete[] table;
  if(dt!=NULL) delete[] dt;
}

// Constructor
powerlaw::powerlaw(double _alpha, int _mini, int _maxi) {
  alpha = _alpha;
  mini = _mini;
  maxi = _maxi;
  if(alpha<=2.0 && maxi<0)
    igraph_warningf("powerlaw exponent %f should be > 2 when no "
		    "Maximum is specified", __FILE__, __LINE__, -1, alpha);
  if(alpha<=1.0 && maxi>=0)
    igraph_warningf("powerlaw exponent %f should be > 1", __FILE__, __LINE__,
		    -1, alpha);
  if(maxi>=0 && mini>maxi)
    igraph_warningf("powerlaw max %d should be greater than min %d",
		    __FILE__, __LINE__, -1, maxi, mini);
  table = new int[POWERLAW_TABLE];
  tabulated = 0;
  dt = NULL;
}

// Sample
int powerlaw::sample() {
  if(proba_big!=0 && test_proba(proba_big)) return int(floor(0.5+big_sample(random_float())));
  int r=my_random();
  // table[] contains integer from MY_RAND_MAX downto 0, in blocks. Search block...
  if(r>(MY_RAND_MAX>>max_dt)) return mini;
  int k=0;
  while(k<max_dt) { r<<=1; r+=random_bit(); k++; };
  int a=0;
  int b;
  while((b=dt[k++])<0 || r<table[b]) { 
    if(b>=0) {
      a=b+1;
      if(a==tabulated-1) break;
      r<<=1;
      r+=random_bit();
    }
  }

  // Now that we found the good block, run a dichotomy on this block [a,b]
  while(a<b) {
    int c = (a+b)/2;
    if(r<table[c]) a=c+1;
    else b=c;
  }
  return mini+a;
}

// Proba
double powerlaw::proba(int k) {
  if(k<mini || (maxi>=0 && k>maxi)) return 0.0;
  if(k>=mini+tabulated)
    return proba_big*(big_inv_sample(double(k)-0.5)-big_inv_sample(double(k)+0.5));
  else {
    double div = table_mul;
    int prev_pos_in_table = k-mini-1;
    if(prev_pos_in_table<0) return (double(MY_RAND_MAX)+1.0-double(table[0]>>max_dt))*div;
    // what block are we in ?
    int k=0;
    while(k<max_dt) { div*=0.5; k++; };
    while(dt[k]<0 || dt[k]<prev_pos_in_table) { k++; div*=0.5; };
    double prob2 = double(table[prev_pos_in_table+1]);
    if(dt[k]==prev_pos_in_table) do prob2*=0.5;while(dt[++k]<0);
    return (double(table[prev_pos_in_table])-prob2)*div;
  }
}

// Relative Error
double powerlaw::error() {
  return 1.0/(double(tabulated)*double(tabulated));
}

// Mean
double powerlaw::mean() {
  double sum = 0.0;
  for(int i=mini+tabulated; --i>=mini; ) sum+=double(i)*proba(i);
  // add proba_big * integral(big_sample(t),t=0..1)
  if(proba_big!=0) sum += proba_big*((pow(_a+_b,_exp+1.0)-pow(_b,_exp+1.0))/(_a*(_exp+1.0)) +double(mini)-offset-sum);
  return sum;
}

// Median. Returns integer Med such that P(X<=Med) >= 1/2
int powerlaw::median() {
  if(proba_big>0.5) return int(floor(0.5+big_sample(1.0-0.5/proba_big)));
  double sum = 0.0;
  int i=mini;
  while(sum<0.5) sum+=proba(i++);
  return i-1;
}  

void powerlaw::init_to_offset(double _offset, int _tabulated) {
  offset = _offset;
  tabulated = _tabulated;
  if(maxi>=0 && tabulated > maxi-mini) tabulated=maxi-mini+1;
  double sum = 0.0;
  double item = double(tabulated)+offset;
  // Compute sum of tabulated probabilities
  for(int i=tabulated; i--; ) sum += pow(item-=1.0, -alpha);
  // Compute others parameters : proba_big, table_mul, _a, _b, _exp
  if(maxi>0 && maxi<=mini+tabulated-1) {
    proba_big = 0;
    table_mul = inv_RANDMAX;
  }
  else {
    if(maxi<0) _b = 0.0;
    else _b = pow(double(maxi-mini)+0.5+offset, 1.0-alpha);
    _a = pow(double(tabulated)-0.5+offset,1.0-alpha) - _b;
    _exp = 1.0 / (1.0 - alpha);
    double sum_big = _a*(-_exp);
    proba_big = sum_big / (sum + sum_big);
    table_mul = inv_RANDMAX * sum / (sum + sum_big);
  }
  // How many delimiters will be necessary for the table ?
  max_dt = max(0,int(floor(alpha*log(double(tabulated))/log(2.0)))-6);
  if(dt!=NULL) delete[] dt;
  dt = new int[max_dt+1];
  // Create table as decreasing integers from MY_RAND_MAX+1 (in virtual position -1) down to 0
  // Every time the index crosses a delimiter, numbers get doubled.
  double ssum = 0;
  double mul = (double(MY_RAND_MAX)+1.0)*pow(2.0,max_dt)/sum;
  item = double(tabulated)+offset;
  int k = max_dt;
  dt[k--]=tabulated-1;
  for(int i=tabulated; --i>0; ) {
    table[i] = int(floor(0.5+ssum));
    ssum += mul * pow(item-=1.0,-alpha);
    if(ssum>double(MY_RAND_MAX/2) && k>=0) {
      while((ssum*=0.5)>double(MY_RAND_MAX/2)) { mul*=0.5; dt[k--]=-1; };
      mul*=0.5; dt[k--]=i-1;
    }
  }
  table[0] = int(floor(0.5+ssum));
  max_dt = k+1;
}

void powerlaw::adjust_offset_mean(double _mean, double err, double factor) {
  // Set two bounds for offset
  double ol = offset;
  double oh = offset;
  if(mean()<_mean) {
    do {
      ol = oh;
      oh *= factor;
      init_to_offset(oh, tabulated);
    } while(mean()<_mean);
  }
  else {
    do {
      oh = ol;
      ol /= factor;
      init_to_offset(ol, tabulated);
    } while(mean()>_mean);
  }
  // Now, dichotomy
  while(fabs(oh-ol) > err*ol) {
    double oc = sqrt(oh*ol);
    init_to_offset(oc, tabulated);
    if(mean()<_mean) ol = oc;
    else oh = oc;
  }
  init_to_offset(sqrt(ol*oh), tabulated);
}

double powerlaw::init_to_mean(double _mean) {
  if(maxi>=0 && _mean >= 0.5*double((mini+maxi))) {
    igraph_errorf("Fatal error in powerlaw::init_to_mean(%f): "
		  "Mean must be in ]min, (min+max)/2[ = ]%d, %d[",
		  __FILE__, __LINE__, IGRAPH_EINVAL, 
		  _mean, mini, (mini+maxi)/2);
    return(-1.0);
  }
  init_to_offset(_mean-double(mini), 100);
  adjust_offset_mean(_mean, 0.01, 2);
  init_to_offset(offset, POWERLAW_TABLE);
  double eps = 1.0/(double(POWERLAW_TABLE));
  adjust_offset_mean(_mean, eps*eps, 1.01);
  return offset;
}

} // namespace gengraph
